{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.03819,
     "end_time": "2021-02-13T18:19:43.324200",
     "exception": false,
     "start_time": "2021-02-13T18:19:43.286010",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "\n",
    "\n",
    "This notebook contains an example for teaching.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_execution_state": "idle",
    "_uuid": "051d70d956493feee0c6d64651c6a088724dca2a",
    "papermill": {
     "duration": 0.036479,
     "end_time": "2021-02-13T18:19:43.396666",
     "exception": false,
     "start_time": "2021-02-13T18:19:43.360187",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# A Simple Case Study using Wage Data from 2015"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.036639,
     "end_time": "2021-02-13T18:19:43.468425",
     "exception": false,
     "start_time": "2021-02-13T18:19:43.431786",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We illustrate how to predict an outcome variable Y in a high-dimensional setting, where the number of covariates $p$ is large in relation to the sample size $n$. So far we have used linear prediction rules, e.g. Lasso regression, for estimation.\n",
    "Now, we also consider nonlinear prediction rules including tree-based methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.034705,
     "end_time": "2021-02-13T18:19:43.537814",
     "exception": false,
     "start_time": "2021-02-13T18:19:43.503109",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.036082,
     "end_time": "2021-02-13T18:19:43.609347",
     "exception": false,
     "start_time": "2021-02-13T18:19:43.573265",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Again, we consider data from the U.S. March Supplement of the Current Population Survey (CPS) in 2015.\n",
    "The preproccessed sample consists of $5150$ never-married individuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.279387,
     "end_time": "2021-02-13T18:19:43.923823",
     "exception": false,
     "start_time": "2021-02-13T18:19:43.644436",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "load(\"../input/wage2015-inference/wage2015_subsample_inference.Rdata\")\n",
    "dim(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.034902,
     "end_time": "2021-02-13T18:19:43.994834",
     "exception": false,
     "start_time": "2021-02-13T18:19:43.959932",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The outcomes $Y_i$'s are hourly (log) wages of never-married workers living in the U.S. The raw regressors $Z_i$'s consist of a variety of characteristics, including experience, education and industry and occupation indicators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.091723,
     "end_time": "2021-02-13T18:19:44.123394",
     "exception": false,
     "start_time": "2021-02-13T18:19:44.031671",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "Z <- subset(data,select=-c(lwage,wage)) # regressors\n",
    "colnames(Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.037074,
     "end_time": "2021-02-13T18:19:44.196749",
     "exception": false,
     "start_time": "2021-02-13T18:19:44.159675",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The following figure shows the weekly wage distribution from the US survey data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.443391,
     "end_time": "2021-02-13T18:19:44.677379",
     "exception": false,
     "start_time": "2021-02-13T18:19:44.233988",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "hist(data$wage, xlab= \"hourly wage\", main=\"Empirical wage distribution from the US survey data\", breaks= 35)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.036602,
     "end_time": "2021-02-13T18:19:44.752465",
     "exception": false,
     "start_time": "2021-02-13T18:19:44.715863",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Wages show a high degree of skewness. Hence, wages are transformed in almost all studies by\n",
    "the logarithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.036009,
     "end_time": "2021-02-13T18:19:44.826260",
     "exception": false,
     "start_time": "2021-02-13T18:19:44.790251",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.036925,
     "end_time": "2021-02-13T18:19:44.899159",
     "exception": false,
     "start_time": "2021-02-13T18:19:44.862234",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Due to the skewness of the data, we are considering log wages which leads to the following regression model\n",
    "\n",
    "$$log(wage) = g(Z) + \\epsilon.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.036183,
     "end_time": "2021-02-13T18:19:44.971528",
     "exception": false,
     "start_time": "2021-02-13T18:19:44.935345",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We will estimate the two sets of prediction rules: Linear and Nonlinear Models.\n",
    "In linear models, we estimate the prediction rule of the form\n",
    "\n",
    "$$\\hat g(Z) = \\hat \\beta'X.$$\n",
    "Again, we generate $X$ in two ways:\n",
    " \n",
    "1. Basic Model:   $X$ consists of a set of raw regressors (e.g. gender, experience, education indicators, regional indicators).\n",
    "\n",
    "\n",
    "2. Flexible Model:  $X$ consists of all raw regressors from the basic model plus occupation and industry indicators, transformations (e.g., ${exp}^2$ and ${exp}^3$) and additional two-way interactions.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.037318,
     "end_time": "2021-02-13T18:19:45.044959",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.007641",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "To evaluate the out-of-sample performance, we split the data first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.062188,
     "end_time": "2021-02-13T18:19:45.143118",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.080930",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "set.seed(1234)\n",
    "training <- sample(nrow(data), nrow(data)*(3/4), replace=FALSE)\n",
    "\n",
    "data_train <- data[training,]\n",
    "data_test <- data[-training,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.038774,
     "end_time": "2021-02-13T18:19:45.217757",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.178983",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We construct the two different model matrices $X_{basic}$ and $X_{flex}$ for both the training and the test sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.094135,
     "end_time": "2021-02-13T18:19:45.347955",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.253820",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "X_basic <-  \"sex + exp1 + exp2+ shs + hsg+ scl + clg + mw + so + we + occ2+ ind2\"\n",
    "X_flex <- \"sex + exp1 + exp2 + shs+hsg+scl+clg+occ2+ind2+mw+so+we + (exp1+exp2+exp3+exp4)*(shs+hsg+scl+clg+occ2+ind2+mw+so+we)\"\n",
    "formula_basic <- as.formula(paste(\"lwage\", \"~\", X_basic))\n",
    "formula_flex <- as.formula(paste(\"lwage\", \"~\", X_flex))\n",
    "\n",
    "model_X_basic_train <- model.matrix(formula_basic,data_train)\n",
    "model_X_basic_test <- model.matrix(formula_basic,data_test)\n",
    "p_basic <- dim(model_X_basic_train)[2]\n",
    "model_X_flex_train <- model.matrix(formula_flex,data_train)\n",
    "model_X_flex_test <- model.matrix(formula_flex,data_test)\n",
    "p_flex <- dim(model_X_flex_train)[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.060969,
     "end_time": "2021-02-13T18:19:45.445389",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.384420",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "Y_train <- data_train$lwage\n",
    "Y_test <- data_test$lwage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.062723,
     "end_time": "2021-02-13T18:19:45.545189",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.482466",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "p_basic\n",
    "p_flex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.037704,
     "end_time": "2021-02-13T18:19:45.622370",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.584666",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "As known from our first lab, the basic model consists of $10$ regressors and the flexible model of $246$ regressors. Let us fit our models to the training sample using the two different model specifications. We are starting by running a simple ols regression. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.038763,
     "end_time": "2021-02-13T18:19:45.699126",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.660363",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### OLS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.039458,
     "end_time": "2021-02-13T18:19:45.779460",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.740002",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We fit the basic model to our training data by running an ols regression and compute the mean squared error on the test sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.069537,
     "end_time": "2021-02-13T18:19:45.887169",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.817632",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# ols (basic model)\n",
    "fit.lm.basic <- lm(formula_basic, data_train)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.074423,
     "end_time": "2021-02-13T18:19:45.999870",
     "exception": false,
     "start_time": "2021-02-13T18:19:45.925447",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Compute the Out-Of-Sample Performance\n",
    "yhat.lm.basic <- predict(fit.lm.basic, newdata=data_test)\n",
    "cat(\"The mean squared error (MSE) using the basic model is equal to\" , mean((Y_test-yhat.lm.basic)^2)) # MSE OLS (basic model)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.052764,
     "end_time": "2021-02-13T18:19:46.122829",
     "exception": false,
     "start_time": "2021-02-13T18:19:46.070065",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "To determine the out-of-sample $MSE$ and the standard error in one step, we can use the function *lm*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.076484,
     "end_time": "2021-02-13T18:19:46.239015",
     "exception": false,
     "start_time": "2021-02-13T18:19:46.162531",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "MSE.lm.basic <- summary(lm((Y_test-yhat.lm.basic)^2~1))$coef[1:2]\n",
    "MSE.lm.basic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.039088,
     "end_time": "2021-02-13T18:19:46.317915",
     "exception": false,
     "start_time": "2021-02-13T18:19:46.278827",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We also compute the out-of-sample $R^2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.057098,
     "end_time": "2021-02-13T18:19:46.413754",
     "exception": false,
     "start_time": "2021-02-13T18:19:46.356656",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "R2.lm.basic <- 1-MSE.lm.basic[1]/var(Y_test)\n",
    "cat(\"The R^2 using the basic model is equal to\",R2.lm.basic) # MSE OLS (basic model) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.039585,
     "end_time": "2021-02-13T18:19:46.492903",
     "exception": false,
     "start_time": "2021-02-13T18:19:46.453318",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We repeat the same procedure for the flexible model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.198636,
     "end_time": "2021-02-13T18:19:46.730717",
     "exception": false,
     "start_time": "2021-02-13T18:19:46.532081",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# ols (flexible model)\n",
    "fit.lm.flex <- lm(formula_flex, data_train)  \n",
    "# Compute the Out-Of-Sample Performance\n",
    "options(warn=-1)\n",
    "yhat.lm.flex <- predict(fit.lm.flex, newdata=data_test)\n",
    "MSE.lm.flex <- summary(lm((Y_test-yhat.lm.flex)^2~1))$coef[1:2]\n",
    "R2.lm.flex <- 1-MSE.lm.flex[1]/var(Y_test)\n",
    "cat(\"The R^2 using the flexible model is equal to\",R2.lm.flex) # MSE OLS (flexible model) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.051953,
     "end_time": "2021-02-13T18:19:46.853182",
     "exception": false,
     "start_time": "2021-02-13T18:19:46.801229",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We observe that ols regression works better for the basic model with smaller $p/n$ ratio. We are proceeding by running lasso regressions and its versions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.042521,
     "end_time": "2021-02-13T18:19:46.935859",
     "exception": false,
     "start_time": "2021-02-13T18:19:46.893338",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Lasso, Ridge and Elastic Net\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.040161,
     "end_time": "2021-02-13T18:19:47.015626",
     "exception": false,
     "start_time": "2021-02-13T18:19:46.975465",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Considering the basic model, we run a lasso/post-lasso regression first and then we compute the measures for the out-of-sample performance. Note that applying the package *hdm* and the function *rlasso* we rely on a theoretical based choice of the penalty level $\\lambda$ in the lasso regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.577781,
     "end_time": "2021-02-13T18:19:47.634269",
     "exception": false,
     "start_time": "2021-02-13T18:19:47.056488",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# lasso and versions\n",
    "library(hdm) \n",
    "fit.rlasso  <- rlasso(formula_basic, data_train, post=FALSE)\n",
    "fit.rlasso.post <- rlasso(formula_basic, data_train, post=TRUE)\n",
    "yhat.rlasso   <- predict(fit.rlasso, newdata=data_test)\n",
    "yhat.rlasso.post   <- predict(fit.rlasso.post, newdata=data_test)\n",
    "\n",
    "MSE.lasso <- summary(lm((Y_test-yhat.rlasso)^2~1))$coef[1:2]\n",
    "MSE.lasso.post <- summary(lm((Y_test-yhat.rlasso.post)^2~1))$coef[1:2]\n",
    "\n",
    "R2.lasso <- 1-MSE.lasso[1]/var(Y_test)\n",
    "R2.lasso.post <- 1-MSE.lasso.post[1]/var(Y_test)\n",
    "cat(\"The R^2 using the basic model is equal to\",R2.lasso,\"for lasso and\",R2.lasso.post,\"for post-lasso\") # R^2 lasso/post-lasso (basic model) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.049543,
     "end_time": "2021-02-13T18:19:47.757271",
     "exception": false,
     "start_time": "2021-02-13T18:19:47.707728",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Now, we repeat the same procedure for the flexible model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 3.430649,
     "end_time": "2021-02-13T18:19:51.229007",
     "exception": false,
     "start_time": "2021-02-13T18:19:47.798358",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fit.rlasso.flex  <- rlasso(formula_flex, data_train, post=FALSE)\n",
    "fit.rlasso.post.flex <- rlasso(formula_flex, data_train, post=TRUE)\n",
    "yhat.rlasso.flex   <- predict(fit.rlasso.flex, newdata=data_test)\n",
    "yhat.rlasso.post.flex   <- predict(fit.rlasso.post.flex, newdata=data_test)\n",
    "\n",
    "MSE.lasso.flex <- summary(lm((Y_test-yhat.rlasso.flex)^2~1))$coef[1:2]\n",
    "MSE.lasso.post.flex <- summary(lm((Y_test-yhat.rlasso.post.flex)^2~1))$coef[1:2]\n",
    "\n",
    "R2.lasso.flex <- 1-MSE.lasso.flex[1]/var(Y_test)\n",
    "R2.lasso.post.flex <- 1-MSE.lasso.post.flex[1]/var(Y_test)\n",
    "cat(\"The R^2 using the flexible model is equal to\",R2.lasso.flex,\"for lasso and\",R2.lasso.post.flex,\"for post-lasso\") # R^2 lasso/post-lasso (flexible model) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.04206,
     "end_time": "2021-02-13T18:19:51.353816",
     "exception": false,
     "start_time": "2021-02-13T18:19:51.311756",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "It is worth to notice that lasso regression works better for the more complex model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.041452,
     "end_time": "2021-02-13T18:19:51.436401",
     "exception": false,
     "start_time": "2021-02-13T18:19:51.394949",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "In contrast to a theoretical based choice of the tuning parameter $\\lambda$ in the lasso regression, we can also use cross-validation to determine the penalty level by applying the package *glmnet* and the function cv.glmnet. In this context, we also run a ridge and a elastic net regression by adjusting the parameter *alpha*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 2.248453,
     "end_time": "2021-02-13T18:19:53.725885",
     "exception": false,
     "start_time": "2021-02-13T18:19:51.477432",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "library(glmnet)\n",
    "fit.lasso.cv   <- cv.glmnet(model_X_basic_train, Y_train, family=\"gaussian\", alpha=1)\n",
    "fit.ridge   <- cv.glmnet(model_X_basic_train, Y_train, family=\"gaussian\", alpha=0)\n",
    "fit.elnet   <- cv.glmnet(model_X_basic_train, Y_train, family=\"gaussian\", alpha=.5)\n",
    "\n",
    "yhat.lasso.cv    <- predict(fit.lasso.cv, newx = model_X_basic_test)\n",
    "yhat.ridge   <- predict(fit.ridge, newx = model_X_basic_test)\n",
    "yhat.elnet   <- predict(fit.elnet, newx = model_X_basic_test)\n",
    "\n",
    "MSE.lasso.cv <- summary(lm((Y_test-yhat.lasso.cv)^2~1))$coef[1:2]\n",
    "MSE.ridge <- summary(lm((Y_test-yhat.ridge)^2~1))$coef[1:2]\n",
    "MSE.elnet <- summary(lm((Y_test-yhat.elnet)^2~1))$coef[1:2]\n",
    "\n",
    "R2.lasso.cv <- 1-MSE.lasso.cv[1]/var(Y_test)\n",
    "R2.ridge <- 1-MSE.ridge[1]/var(Y_test)\n",
    "R2.elnet <- 1-MSE.elnet[1]/var(Y_test)\n",
    "\n",
    "# R^2 using cross-validation (basic model) \n",
    "cat(\"R^2 using cross-validation for lasso, ridge and elastic net in the basic model:\",R2.lasso.cv,R2.ridge,R2.elnet)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.042613,
     "end_time": "2021-02-13T18:19:53.812553",
     "exception": false,
     "start_time": "2021-02-13T18:19:53.769940",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Note that the following calculations for the flexible model need some computation time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 13.588391,
     "end_time": "2021-02-13T18:20:07.443188",
     "exception": false,
     "start_time": "2021-02-13T18:19:53.854797",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fit.lasso.cv.flex   <- cv.glmnet(model_X_flex_train, Y_train, family=\"gaussian\", alpha=1)\n",
    "fit.ridge.flex   <- cv.glmnet(model_X_flex_train, Y_train, family=\"gaussian\", alpha=0)\n",
    "fit.elnet.flex   <- cv.glmnet(model_X_flex_train, Y_train, family=\"gaussian\", alpha=.5)\n",
    "\n",
    "yhat.lasso.cv.flex    <- predict(fit.lasso.cv.flex , newx = model_X_flex_test)\n",
    "yhat.ridge.flex    <- predict(fit.ridge.flex , newx = model_X_flex_test)\n",
    "yhat.elnet.flex    <- predict(fit.elnet.flex , newx = model_X_flex_test)\n",
    "\n",
    "MSE.lasso.cv.flex  <- summary(lm((Y_test-yhat.lasso.cv.flex )^2~1))$coef[1:2]\n",
    "MSE.ridge.flex  <- summary(lm((Y_test-yhat.ridge.flex )^2~1))$coef[1:2]\n",
    "MSE.elnet.flex  <- summary(lm((Y_test-yhat.elnet.flex )^2~1))$coef[1:2]\n",
    "\n",
    "R2.lasso.cv.flex  <- 1-MSE.lasso.cv.flex [1]/var(Y_test)\n",
    "R2.ridge.flex  <- 1-MSE.ridge.flex [1]/var(Y_test)\n",
    "R2.elnet.flex  <- 1-MSE.elnet.flex [1]/var(Y_test)\n",
    "\n",
    "# R^2 using cross-validation (flexible model) \n",
    "cat(\"R^2 using cross-validation for lasso, ridge and elastic net in the flexible model:\",R2.lasso.cv.flex,R2.ridge.flex,R2.elnet.flex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.04263,
     "end_time": "2021-02-13T18:20:07.529566",
     "exception": false,
     "start_time": "2021-02-13T18:20:07.486936",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The performance of the lasso regression with cross-validated penalty is quite similar to the performance of lasso using a theoretical based choice of the tuning parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.042859,
     "end_time": "2021-02-13T18:20:07.614751",
     "exception": false,
     "start_time": "2021-02-13T18:20:07.571892",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Non-linear models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.042125,
     "end_time": "2021-02-13T18:20:07.699092",
     "exception": false,
     "start_time": "2021-02-13T18:20:07.656967",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Besides linear regression models, we consider nonlinear regression models to build a predictive model. We are applying regression trees, random forests, boosted trees and neural nets to estimate the regression function $g(X)$. First, we load the relevant libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.345105,
     "end_time": "2021-02-13T18:20:08.087641",
     "exception": false,
     "start_time": "2021-02-13T18:20:07.742536",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "library(randomForest)\n",
    "library(rpart)\n",
    "library(nnet)\n",
    "library(gbm)\n",
    "library(rpart.plot)\n",
    "library(keras)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.043418,
     "end_time": "2021-02-13T18:20:08.174622",
     "exception": false,
     "start_time": "2021-02-13T18:20:08.131204",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "and we illustrate the application of regression trees."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.043267,
     "end_time": "2021-02-13T18:20:08.261600",
     "exception": false,
     "start_time": "2021-02-13T18:20:08.218333",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Regression Trees"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.043445,
     "end_time": "2021-02-13T18:20:08.348402",
     "exception": false,
     "start_time": "2021-02-13T18:20:08.304957",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We fit a regression tree to the training data using the basic model. The variable *cp* controls the complexity of the regression tree, i.e. how deep we build the tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 1.843972,
     "end_time": "2021-02-13T18:20:10.235503",
     "exception": false,
     "start_time": "2021-02-13T18:20:08.391531",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# tree\n",
    "fit.trees   <- rpart(formula_basic, data_train,cp = 0.001)\n",
    "prp(fit.trees,leaf.round=1, space=2, yspace=2,split.space=2,shadow.col = \"gray\",trace = 1) # plotting the tree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.046456,
     "end_time": "2021-02-13T18:20:10.328795",
     "exception": false,
     "start_time": "2021-02-13T18:20:10.282339",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "An important method to improve predictive performance is called \"Pruning the Tree\". This\n",
    "means the process of cutting down the branches of a tree. We apply pruning to the complex tree above to reduce the depth. Initially, we determine the optimal complexity of the regression tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.070106,
     "end_time": "2021-02-13T18:20:10.445828",
     "exception": false,
     "start_time": "2021-02-13T18:20:10.375722",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "bestcp<- fit.trees$cptable[which.min(fit.trees$cptable[,\"xerror\"]),\"CP\"]\n",
    "bestcp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.047157,
     "end_time": "2021-02-13T18:20:10.540327",
     "exception": false,
     "start_time": "2021-02-13T18:20:10.493170",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Now, we can prune the tree and visualize the prediction rule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.543483,
     "end_time": "2021-02-13T18:20:11.131455",
     "exception": false,
     "start_time": "2021-02-13T18:20:10.587972",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "fit.prunedtree <- prune(fit.trees,cp=bestcp)\n",
    "prp(fit.prunedtree,leaf.round=1, space=3, yspace=3, split.space=7, shadow.col = \"gray\",trace = 1,yesno=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.052169,
     "end_time": "2021-02-13T18:20:11.234467",
     "exception": false,
     "start_time": "2021-02-13T18:20:11.182298",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "E.g., in the pruned tree the predicted hourly log wage for high-school graduates with more than $9.5$ years of experience is $2.8$, and otherwise is $2.6$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.04994,
     "end_time": "2021-02-13T18:20:11.334448",
     "exception": false,
     "start_time": "2021-02-13T18:20:11.284508",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Finally, we calculate the mean-squared error and the $R^2$ on the test sample to evaluate the out-of-sample performance of the pruned tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.079534,
     "end_time": "2021-02-13T18:20:11.463701",
     "exception": false,
     "start_time": "2021-02-13T18:20:11.384167",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "yhat.pt <- predict(fit.prunedtree,newdata=data_test)\n",
    "MSE.pt <- summary(lm((Y_test-yhat.pt)^2~1))$coef[1:2]\n",
    "R2.pt  <- 1-MSE.pt[1]/var(Y_test)\n",
    "\n",
    "# R^2 of the pruned tree\n",
    "cat(\"R^2 of the pruned tree:\",R2.pt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.052287,
     "end_time": "2021-02-13T18:20:11.566330",
     "exception": false,
     "start_time": "2021-02-13T18:20:11.514043",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Random Forest and Boosted Trees"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.050794,
     "end_time": "2021-02-13T18:20:11.667980",
     "exception": false,
     "start_time": "2021-02-13T18:20:11.617186",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "In the next step, we apply the more advanced tree-based methods random forest and boosted trees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 56.677891,
     "end_time": "2021-02-13T18:21:08.396363",
     "exception": false,
     "start_time": "2021-02-13T18:20:11.718472",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "## Applying the methods\n",
    "# random forest\n",
    "fit.rf       <- randomForest(formula_basic, ntree=2000, nodesize=5, data=data_train)\n",
    "# for tuning: adjust input \"mtry\" to change the number of variables randomly sampled as candidates at each split\n",
    "\n",
    "# boosting\n",
    "fit.boost   <- gbm(formula_basic, data=data_train, distribution= \"gaussian\", bag.fraction = .5, interaction.depth=2, n.trees=1000, shrinkage=.01)\n",
    "best.boost  <- gbm.perf(fit.boost, plot.it = FALSE) # cross-validation to determine when to stop\n",
    "\n",
    "## Evaluating the methods\n",
    "yhat.rf       <- predict(fit.rf, newdata=data_test) # prediction\n",
    "yhat.boost    <- predict(fit.boost, newdata=data_test, n.trees=best.boost)\n",
    "\n",
    "MSE.rf       = summary(lm((Y_test-yhat.rf)^2~1))$coef[1:2]\n",
    "MSE.boost    = summary(lm((Y_test-yhat.boost)^2~1))$coef[1:2]\n",
    "R2.rf  <- 1-MSE.rf[1]/var(Y_test)\n",
    "R2.boost  <- 1-MSE.boost[1]/var(Y_test)\n",
    "\n",
    "# printing R^2\n",
    "cat(\"R^2 of the random forest and boosted trees:\",R2.rf,R2.boost)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.051225,
     "end_time": "2021-02-13T18:21:08.500313",
     "exception": false,
     "start_time": "2021-02-13T18:21:08.449088",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "To conclude, let us have a look at our results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.052403,
     "end_time": "2021-02-13T18:21:08.603976",
     "exception": false,
     "start_time": "2021-02-13T18:21:08.551573",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.167847,
     "end_time": "2021-02-13T18:21:08.823485",
     "exception": false,
     "start_time": "2021-02-13T18:21:08.655638",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "library(xtable)\n",
    "table<- matrix(0, 15, 3)\n",
    "table[1,1:2]   <- MSE.lm.basic\n",
    "table[2,1:2]   <- MSE.lm.flex\n",
    "table[3,1:2]   <- MSE.lasso\n",
    "table[4,1:2]   <- MSE.lasso.post\n",
    "table[5,1:2]   <- MSE.lasso.flex\n",
    "table[6,1:2]   <- MSE.lasso.post.flex\n",
    "table[7,1:2]   <- MSE.lasso.cv\n",
    "table[8,1:2]   <- MSE.ridge\n",
    "table[9,1:2]   <- MSE.elnet\n",
    "table[10,1:2]   <- MSE.lasso.cv.flex\n",
    "table[11,1:2]  <- MSE.ridge.flex\n",
    "table[12,1:2]  <- MSE.elnet.flex\n",
    "table[13,1:2]  <- MSE.rf\n",
    "table[14,1:2]  <- MSE.boost\n",
    "table[15,1:2]  <- MSE.pt\n",
    "\n",
    "\n",
    "\n",
    "table[1,3]   <- R2.lm.basic\n",
    "table[2,3]   <- R2.lm.flex\n",
    "table[3,3]   <- R2.lasso\n",
    "table[4,3]   <- R2.lasso.post\n",
    "table[5,3]   <- R2.lasso.flex\n",
    "table[6,3]   <- R2.lasso.post.flex\n",
    "table[7,3]   <- R2.lasso.cv\n",
    "table[8,3]   <- R2.ridge\n",
    "table[9,3]   <- R2.elnet\n",
    "table[10,3]   <- R2.lasso.cv.flex\n",
    "table[11,3]  <- R2.ridge.flex\n",
    "table[12,3]  <- R2.elnet.flex\n",
    "table[13,3]  <- R2.rf\n",
    "table[14,3]  <- R2.boost\n",
    "table[15,3]  <- R2.pt\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "colnames(table)<- c(\"MSE\", \"S.E. for MSE\", \"R-squared\")\n",
    "rownames(table)<- c(\"Least Squares (basic)\",\"Least Squares (flexible)\", \"Lasso\", \"Post-Lasso\",\"Lasso (flexible)\",\"Post-Lasso (flexible)\", \n",
    "                    \"Cross-Validated lasso\", \"Cross-Validated ridge\",\"Cross-Validated elnet\",\"Cross-Validated lasso (flexible)\",\"Cross-Validated ridge (flexible)\",\"Cross-Validated elnet (flexible)\",  \n",
    "                    \"Random Forest\",\"Boosted Trees\", \"Pruned Tree\")\n",
    "tab <- xtable(table, digits =3)\n",
    "print(tab,type=\"latex\") # set type=\"latex\" for printing table in LaTeX\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.052897,
     "end_time": "2021-02-13T18:21:08.930888",
     "exception": false,
     "start_time": "2021-02-13T18:21:08.877991",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Above, we displayed the results for a single split of data into the training and testing part. The table shows the test MSE in column 1 as well as the standard error in column 2 and the test $R^2$\n",
    "in column 3. We see that the prediction rule produced by Elastic Net using the flexible model performs the best here, giving the lowest test MSE. Cross-Validated Lasso and Ridge, perform nearly as well. For any two of these methods, their testing MSEs are within one standard error of each other. Remarkably, OLS on a simple model performs extremely well, almost as well as best tree based method Random Forest. On the other hand, OLS on a flexible model with many regressors performs very poorly giving the highest test MSE. It is worth to notice that the nonlinear models, e.g. Random Forest, are not tuned. Thus, there is a lot of potential to improve the performance of the nonlinear methods we used in the analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.052594,
     "end_time": "2021-02-13T18:21:09.036009",
     "exception": false,
     "start_time": "2021-02-13T18:21:08.983415",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### Ensemble learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.053134,
     "end_time": "2021-02-13T18:21:09.146558",
     "exception": false,
     "start_time": "2021-02-13T18:21:09.093424",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "In the final step, we can build a prediction model by combing the strengths of the models we considered so far. This ensemble method is of the form\n",
    "\t$$ f(x) = \\sum_{k=1}^K \\alpha_k f_k(x) $$\n",
    "where the $f_k$'s denote our prediction rules from the table above and the $\\alpha_k$'s are the corresponding weights."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.053696,
     "end_time": "2021-02-13T18:21:09.254342",
     "exception": false,
     "start_time": "2021-02-13T18:21:09.200646",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We focus on the prediction rules based on OLS, Post-Lasso, Elastic Net, Pruned Tree, Random Forest, Boosted Trees, and Neural Network and combine these methods into an ensemble method. The weights can be determined by a simple ols regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.079851,
     "end_time": "2021-02-13T18:21:09.388686",
     "exception": false,
     "start_time": "2021-02-13T18:21:09.308835",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "ensemble.ols <- summary(lm(Y_test~ yhat.lm.basic + yhat.rlasso.post.flex + yhat.elnet.flex+ yhat.pt+ yhat.rf + yhat.boost))\n",
    "ensemble.ols"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.054822,
     "end_time": "2021-02-13T18:21:09.498067",
     "exception": false,
     "start_time": "2021-02-13T18:21:09.443245",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Alternatively, we can determine the weights via lasso regression. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.175196,
     "end_time": "2021-02-13T18:21:09.727077",
     "exception": false,
     "start_time": "2021-02-13T18:21:09.551881",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "ensemble.lasso <- summary(rlasso(Y_test~ yhat.lm.basic + yhat.rlasso.post.flex + yhat.elnet.flex+ yhat.pt+ yhat.rf + yhat.boost))\n",
    "ensemble.lasso"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.055874,
     "end_time": "2021-02-13T18:21:09.838636",
     "exception": false,
     "start_time": "2021-02-13T18:21:09.782762",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The estimated weights are shown in the following table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "papermill": {
     "duration": 0.094431,
     "end_time": "2021-02-13T18:21:09.988946",
     "exception": false,
     "start_time": "2021-02-13T18:21:09.894515",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "table<- matrix(0, 7, 2)\n",
    "table[1:7,1]   <- ensemble.ols$coef[1:7]\n",
    "table[1:7,2]   <- ensemble.lasso$coef[1:7]\n",
    "\n",
    "\n",
    "colnames(table)<- c(\"Weight OLS\", \"Weight Lasso\")\n",
    "rownames(table)<- c(\"Constant\",\"Least Squares (basic)\",\"Post-Lasso (flexible)\", \"Cross-Validated elnet (flexible)\", \"Pruned Tree\",\n",
    "                    \"Random Forest\",\"Boosted Trees\")\n",
    "tab <- xtable(table, digits =3)\n",
    "print(tab,type=\"latex\") # set type=\"latex\" for printing table in LaTeX\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "papermill": {
     "duration": 0.056002,
     "end_time": "2021-02-13T18:21:10.101284",
     "exception": false,
     "start_time": "2021-02-13T18:21:10.045282",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Further, the $R^2$ for the test sample gets improved from $30\\%$ obtained by OLS to about $31\\%$ obtained by the ensemble method. We see that it is very powerful to aggregate prediction rules into an ensemble rule. Nevertheless, it is worth to notice that we should compare the ensemble method and the single rules on an additional validation set to ensure a fair comparison."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 90.376935,
   "end_time": "2021-02-13T18:21:10.266455",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2021-02-13T18:19:39.889520",
   "version": "2.2.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
